number.of.neutral.patterns <- NULL
sw <- 0
# Sweep div prob
for (divprob in 0:4){
    # Sweep selective coefficient
    for (selcoeff in 0:5){
	# Sweep selective mut rate
	for (selmu in 0:3){
	    # Sweep neutral mut rate
	    for (neutmu in 0:2){
		# Sweep replicates
		for (replicates in 0:2){
		    sw <- sw + 1
		    # Load state of grid at end of simulation
		    grid.log <- read.csv(paste("../sim",sw,"/sim",sw,".grid.log",sep=""),header=FALSE)		    
		    grid.log <- grid.log[grid.log[,1]==7300,]
		    # -1 for NA's corresponding to dead crypts
		    number.of.neutral.patterns[sw] <- length(unique(grid.log[,5])-1)
		    print (paste("Sim",sw,"done."))
		}
	    }
	}	
    }
}

write.table(number.of.neutral.patterns,"simulation.number.and.number.of.neutral.patterns.csv",quote=FALSE,sep=",",col.names=FALSE,row.names=FALSE)


library(gplots)

# Read no of neut pat per sim
number.of.neutral.patterns <- read.csv("simulation.number.and.number.of.neutral.patterns.csv",header=FALSE)

divproblab <- c("r=0","r=0.25","r=0.5","r=0.75","r=1")
selcoefflab <- c("s=0.001","s=0.01","s=0.1","s=0.5","s=1","s=2")
selmulab <- c("smu=10^-5","smu=10^-6","smu=10^-7","smu=10^-8")
neutmulab <- c("nmu=10^-4","nmu=10^-5","nmu=10^-7")

res <- matrix(nrow=72,ncol=5)
sw <- 0
for (divprob in 1:5){
selcoeffrow <- 1
for (selcoeff in 1:6){
    for (selmu in 1:4){    	
    	neutmurow <- 1
    	for (neutmu in 1:3){
    	    sw <- sw + 1
res[selcoeffrow-1+neutmurow-1+selmu,divprob] <- mean(c(number.of.neutral.patterns[sw,1],number.of.neutral.patterns[sw+1,1],number.of.neutral.patterns[sw+2,1]))
    	    sw <- sw + 2	    
       	    neutmurow <- neutmurow + 4
    	}
    }
    selcoeffrow <- selcoeffrow + 12
}
}

res <- round(res)


pdf("number.of.neutral.patterns.pdf",width=6, height=6)
sw <- 0
for (selcoeff in 1:6){
for (neutmu in 1:3){
sw <- sw + 1
subres <- res[(sw*4-3):(sw*4),] 
colnames(subres) <- divproblab
rownames(subres) <- selmulab

mp <- barplot2(subres, beside = TRUE,
        col = gray(1-(0:3/3)),        
        main = paste(selcoefflab[selcoeff],neutmulab[neutmu]), font.main = 4,
        plot.grid = TRUE)
#legend(0,max(subres)*9/10,legend=rownames(subres),fill=gray(1-(0:5/5)))
mtext(side=2, text="Average number of patterns per 3 runs",line=2)
mtext(side=1, text="Crypt division probability (r)",line=2)
}
}
dev.off()

